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The various approximations used in the construction of a first-principles effective Hamiltonian 
for BaTiOa, and their effects on the calculated transition temperatures, are discussed. An effective 
Ifamiltonian for BaTiOa is constructed not from first-principles calculations, but from the struc- 
tural energetics of an atomistic shell model for BaTiOa of Tinte et al. This allows the elimination of 
certain uncontrolled approximations that arise in the comparison of first-principles effective Hamil- 
tonian results with experimental values and the quantification of errors associated with the selection 
of the effective Hamiltonian subspace and subsequent projection. The discrepancies in transition 
temperatures computed in classical simulations for this effective Hamiltonian and for the atomistic 
shell model are shown to be associated primarily with a poor description of the thermal expansion 
in the former case. This leads to specific proposals for refinements to the first-principles effective 
Hamiltonian method. Our results suggest that there are at least two significant sources of error 
in the effective-Hamiltonian treatment of BaTiOa in the literature, i.e., the improper treatment of 
thermal expansion, and the errors inherent in the first-principles approach itself. 
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I. INTRODUCTION 

First-principles methods constitu.te a powerful tool for 
the study of ferroelectric systems.El Ground state struc- 
tures, phonons, spontaneous polarization, and related 
properties, including piezoelectric and dielectric tensors, 
\ have been accurately calculated for a wide variety of per- 
ovskite oxides as well as other ferroelectric compounds. 

Despite advances in algorithms and computer hard- 
ware, the direct calculation of finite temperature be- 
havior, particularly phase transitions, is still far be- 
yond reach, as such calculations involve thousands of 
atoms. However, indirect methods have been developed 
and apnlied to a lacge number oLgystems, including 
BaTiOa, M PbTi03,H and KN.b03,m and even solid so- 
: lutions like Pb(ZipTii_^)03,Eiy Pb(Sco.5Nbo.5)03,N and 
■ K(Nb^Tai_^)03.El In Refs. g % and interatomic 
\ "shell-model" potentials were parametrized by fitting to 
first-principles results, and finite-temperature behavior 
studied by direct simulation of atomistic systems with 
forces and energies obtained from these potentials. The 
; results in Refs. |,|,| and were obtained by an 

effective Hamiltonian construction in which the full sys- 
tem is mapped by a subspace projection onto a statisti- 
cal mechanical model, with parameters determined from 
first- principles calculations of total energies for small dis- 
tortions of an ideal crystal with the cubic perovskite 
structure. The simple form of the resulting effective 
Hamiltonian allows very-large-scale simulations and aids 
in the conceptual interpretation of the results. The two 
approaches have achieved comparable success in repro- 
ducing many essential features of the phase,|transitions 
of ferroelectrics. For example, for BaTi03,crBl3 the ex- 
perimentally observed cubic-tetragonal-orthorhombic- 
rhombohedral phase sequence is correctly reproduced. 
However, while the experimental values of the transition 



temperatures are 403, 278, and 183 K, classical Monte- 
Carlo simulations of the first-principles effective Hamil- 
tonian give 300, 230, and 200 K, while the shell-model 
results are 210, 135 and lOOK, respectively. In both cases 
a correction for the local-density approximation (LDA) 
lattice constant underestimate was included in the model. 

Understanding the origin of this discrepancy with ex- 
periment may help in the development of improved the- 
oretical methods for the calculation of finite tempera- 
ture behavior. Here, we focus our attention on the first- 
principles effective Hamiltonian approach. The discrep- 
ancies in the transition temperatures could result from 
separate errors introduced at various steps of the analy- 
sis. We correspondingly classify the errors into five types. 
Errors in the configuration energies obtained from first- 
principles calculations will be designated Type I. These 
generally can be systematically reduced, with the excep- 
tion of the uncontrolled approximation in the exchange- 
correlation functional required for the practical imple- 
mentation of density functional theory ("LDA error"). 
Type II errors result from the identification of the rele- 
vant degrees of freedom and the projection and approx- 
imate representation of the effective Hamiltonian in the 
corresponding subspace, and will be the main focus of our 
investigation. Errors in the statistical-mechanical simu- 
lations (Type HI) include finite-size effects and sampling 
errors. In the effective Hamiltonian studies to date, it 
has been feasible to make these errors relatively negligi- 
ble. The importance of Type IV errors, resulting from 
the classical treatment of the ions neglecting quantum 
fiuctuations, has been highlighted in a recent study of 
BaTi03 by Iniguez and Vanderbilt .tj The results of this 
study indicate that the classical approximation raises the 
transition temperatures. Thus, this is not the origin of 
the underestimate for BaTi03. In fact, a correct, fully 
quantum-mechanical treatment would increase, not de- 
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crease, the transition-temperature discrepancy. Finally, 
we note that the experimental samples, even in thermo- 
dynamic equilibrium, contain defects and local nonstoi- 
chiometry which lead to deviations of observed properties 
from those of the assumed ideal crystals (Type V errors). 
These crystal imperfections can have various effects on 
the transition temperatures that are in general difficult 
to model. 

To separate and quantify the role of the various er- 
rors in producing the observed discrepancies in transi- 
tion temperatures, several different approaches could be 
applied. The analysis of Type III and Type IV errors 
has been discussed in the previous paragraph. One way 
to investigate Type I errors would be completely to redo 
the effective Hamiltonian study replacing the LDA with 
a generalized-gradient approximation (GGA). However, 
the latter has not been found to give systematic improve- 
ment in the overall agreement of calculated properties 
with experiment ,E3 and thus the value of such a labor- 
intensive investigation is unclear. In principle. Type II 
errors could be eliminated by comparing the effective- 
Hamiltonian transition temperatures with those obtained 
in a fully ab-initio molecular dynamics or Monte Carlo 
calculation. However, as noted above, doing this type 
of direct calculation for sufficiently large systems is so 
computationally demanding that it is impossible in prac- 
tice even for benchmarking purposes, and calculations for 
small supercells and with reduced sampling would intro- 
duce significant finite-size and statistical errors. 

In this paper we develop and carry out an alternative 
method of isolating and quantifying Type II errors, al- 
lowing us to discuss possible refinements of the effective 
Hamiltonian method to reduce or eliminate them. We 
use the total energies computed with the BaTiOs "shell 
model" interatomic potential of Tinte et ala to con- 
struct an effective Hamiltonian, and compare the com- 
puted transition temperatures with those obtained in di- 
rect classical simulations for the "shell model" system. 
In this comparison, we completely eliminate errors of 
Types I, IV and V, and can easily make errors of Type 
III negligible. Thus, we can attribute any discrepancies 
directly to errors of Type II. While such errors will not be 
quantitatively identical to the corresponding errors made 
in the construction of the effective Hamiltonian directly 
from ab-initio results, the general accuracy and physical 
faithfulness of the shell model interatomic potential to 
BaTiOa should render conclusions based on this analysis 
quite meaningful. 

The paper is organized as follows. Section II pro- 
vides technical details of the BaTiOa shell-model inter- 
atomic potential of Tinte et al. that serves as our ref- 
erence system. In Sec. HI we describe the construc- 
tion of the effective Hamiltonian and the parameters 
obtained by fitting to the shell model, paying special 
attention to the approximations and technicalities in- 
volved. In Section IV we present the results obtained 
from the effective-Hamiltonian and shell-model classi- 
cal statistical-mechanical simulations. The discrepancies 



are analyzed in Section V, and possible improvements 
on the various effective-Hamiltonian approximations are 
discussed. Section VI is devoted to the specific issue 
of modeling the thermal expansion within the effective- 
Hamiltonian approach. Finally, in Section VII we present 
a discussion of the broader implications of our analysis; 
in particular, we speculate on the relative importance 
of errors of Types I, II, and IV in the first-principles 
effective-Hamiltonian treatments currently in the litera- 
ture. 



II. SHELL-MODEL INTERATOMIC 
POTENTIAL 

Of the various types of interatomic potentials, shell 
models are uniquely well suited to giving a good de- 
scription of the lattice dynamics of perovskite oxides. 
The form of the shell-model potential developed for 
BaTiOa in Ref. ^ incorporates earlier observations that 
the oxygen ahgll -core interaction should be nonlinear and 
anisotropic. tllta Each ion (Ba, Ti or O) is modeled as a 
massive core linked to a massless shell. The core-shell 
interactions for Ba and Ti are harmonic and isotropic. 
An anisotropic core-shell interaction is considered at the 
O^^ ions, with a fourth-order core-shell interaction along 
the 0-Ti bond. In addition to the Coulomb interactions 
between ion cores and shells, the model contains pair- 
wise short-range inter-shell potentials of the Buckingham 
type, i.e., V{r) = a exp{—r/p) — c/r^. The Born-Mayer 
form (c = 0) is sufficient for the Ti-0 and Ba-0 short- 
range potential, while for the 0-0 potential the value of 
c is nonzero. The physically important nonlinearities of 
the interatomic interactions are thus naturally incorpo- 
rated into the form of the potential. 

The material-specific parameters in the interatomic 
potential were determined by adjusting them to fit se- 
lected first-principles results computed using the lin- 
earized augmented planewave (LAPW) method. It 
should be noted, however, that the equilibrium lattice 
constant of the cubic phase is fitted to the experimental 
cubic-phase lattice constant extrapolated to K (3.995 
A), not the LAPW lattice constant. The double wells 
for polar distortions along (001), (Oil), and (111) are 
satisfactorily reproduced, as are the phonon dispersion 
curves for the cubic structure at the experimental lat- 
tice constant. The bulk modulus of the cubic phase and 
the anomalous Born effective charges are also in good 
agreement with the first-principles results. Reference ^ 
contains further details about the construction of the in- 
teratomic potential and values of the parameters. 

Finite-temperature properties of the system de- 
scribed by this interatomic potential are investigated by 
constant-pressure molecular dmamics (MD) simulations 
using the DL-POLY package,EEl where the adiabatic dy- 
namics of the electronic shells are approximated by as- 
signing small masses to them. A Hoover constant-((7,T) 
algorithm with external stress set to zero is employed; all 
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cell lengths and cell angles are allowed to fluctuate. The 
time step is 0.4 fs and the total time of each simulation, 
after 2 ps of thermalization, is 20 ps. Results for a 7 x 7 x 7 
periodic supercell (1715 ions plus 1715 shells which are 
additional degrees of freedom) were reported in Rcf. ^. 
It was shown that the cubic-tetragonal-orthorhombic- 
rhombohedral phase sequence is correctly reproduced. 
Good agreement with experimental data was obtained for 
the structural parameters in the various phases as well as 
the volume thermal expansion coefficient, showing that 
the most important nonlinearities have been included in 
the model. However, the transition temperatures are 
rather low compared to experiment (190, 120, and 90 K). 
This discrepancy does not affect the present analysis of 
Type II errors. In the present work, we have expanded 
the supercell to 10x10x10 primitive cells (10000 degrees 
of freedom). This yields essentially the same results, ex- 
cept that the calculated transition temperatures increase 
shghtly (210, 135, and 100 K). Additional MD simula- 
tions were performed at constant volume, using a modi- 
fied Hoover constant- ((t,T) algorithm that allows for fluc- 
tuations in the cell shape. 



III. CONSTRUCTION OF THE EFFECTIVE 
HAMILTONIAN 

In this Section we describe the effective Hamilto- 
nian that we have constructed using the shell model for 
BaTiOa of Tinte et alJa as our target system. The form 
of the effective Hamiltonian is identical to that proposed 
by Zhong et al.p except that the inhomogeneous strain 
variables found to be unimportant in that study are not 
included here. 

An effective Hamiltonian is a Taylor expansion of the 
energy surface of the system around a high-symmetry 
phase in terms of a set of relevant degrees of freedom. For 
ferroelectric perovskites, the most convenient reference 
structure is the cubic paraelectric phase. The relevant de- 
grees of freedom can be identified by studying the energy 
changes induced by small (harmonic) perturbations of 
the reference structure. The low-energy, typically unsta- 
ble, distortions are the relevant ones, and are expressed, 
in the form of local modes or lattice Wannier functions.tZl 
The relevant local modes are those that add up to pro- 
duce the distorted ferroelectric ground-state structure. 
Also, we take the homogeneous strains as relevant and 
include them in the Hamiltonian. 

There are two possible ways of performing the har- 
monic analysis that leads to the identification and cal- 
culation of the relevant lattice Wannier functions. One 
can study either the force-constant matrix (the matrix of 
second derivatives of the energy with respect to atomic 
displacements) or the corresponding dynamical matrix. 
While the former choice leads to a better description of 
the lowest-energy states of the system, the latter pro- 
vides a kinetic decoupling between the relevant and ir- 
relevant (i.e., not considered in the Hamiltonian) degrees 



TABLE I: Expansion parameters of the effective Hamiltonian 
fitted to the shell-model BaTiOa target system. The notation 
is taken from Ref. H. AH the parameters are in atomic units. 



On-site 


AC2 


0.0562 


a 


0.805 


7 


-0.849 




jl 


-0.01424 


h 


-0.01506 






Intersite 


js 


0.00422 




-0.00240 


js 


0.01956 




J6 


0.00100 


h 


0.00050 






Elastic 


Bu 


5.42 


Bi2 


2.06 


B44 


2.07 


Coupling 




-4.16 


B\yy 


-1.19 




-0.44 


Dipole 


z* 


8.153 




5.24 





of freedom. Here we have worked with the force-constant 
matrix, which is more appropriate for the study of equi- 
librium properties. In any case, we find numerically that 
the force-constant and dynamical matrix descriptions are 
essentially identical. ■— ■ 

Once a relevant set of phonon brancheal3 has been 
identified, the calculation of the corresponding lattice 
Wannier functions can be done at different levels of ap- 
proximation. At the crudest level, they can be con- 
structed from phonons at a single k point ,0 more sophisti- 
cated schemes allowing for better descriptions, of the rel- 
evant phonons throughout the Brillouin zonc.EZHlj Here, 
we calculate the local modes from the unstable phonons 
at r (these generate the ferroelectric structure) choosing 
the local modes to be centered on the Ti atom. The re- 
sulting local modes reproduce the unstable phonons at 
zone-boundary points M and R with an accuracy above 
97%. We thus do not expect a better local-mode defini- 
tion to constitute a significant improvement of our effec- 
tive Hamiltonian of shell- model BaTiOs. 

Let Uia denote the local mode amplitude in unit cell 
i along Cartesian direction a. Let rji denote the strains, 
where I is a Voigt index. Our effective Hamiltonian can 
then be written as 

Following Ref. ^, we have written the effective Hamilto- 
nian as the sum of four terms: the on-site self-energy 
of the local modes i?*"^'^, the long-range dipole-dipole 
interactions between local modes E'^^^, the short-range 
interactions between local modes E^^°'-'^, the elastic en- 
ergy E°^^^, and the interaction between local modes and 
strains E''"* . The dependence of each term on the model 
variables is indicated in Eq. (p. It should be noted that 
the form of the Hamiltonian is greatly simplified by the 
cubic symmetry of the reference structure. For example, 
there are no odd terms in Uia- 

The relevant phonon branches are described by the 
harmonic terms in E^^^^ , i?^P', and E^^°^'^. Anharmonic 
terms for the local modes, required to stabilize the low- 
symmetry phases, are included only in E^^^^ (the "local- 
anharmonicity approximation"). In both E'^^'^ and i?"^* 
only the lowest-order terms in the expansion are consid- 
ered. This constitutes the minimal microscopic model for 
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FIG. 1: Behavior of three Cartesian components of the mean 
(super-cell averaged) polarization of BaTiOa as obtained from 
(a) MC simulations using the effective Hamiltonian, and (b) 
MD simulations directly from the shell model. 

the description of ferroelectricity in BaTiOs. Including 
higher-order terms in the Hamiltonian would constitute 
a systematic improvement of the model. 

As already mentioned, the Hamiltonian we have just 
described is essentially that of Ref. ^, except the param- 
eters were fitted to a series of shell-model calculations 
of total energies and forces, instead of ab-initio results, 
for a set of distorted structures. Following the notation 
of Ref. ^, we list in Tabic || the parameters defining our 
effective Hamiltonian for shell- model BaTiOa. 

Monte Carlo (MC) simulations were performed to cal- 
culate the finite-temperature properties of the Hamilto- 
nian. We simulated a 10x10x10 supercell with periodic 
boundary conditions, and typically did 30,000 MC sweeps 
to equilibrate the system and 50,000 sweeps more to ob- 
tain averages of local-mode variables with a statistical er- 
ror below 10%. The temperature was increased in small 
steps of 5 K. We monitored the behavior of the homo- 
geneous strain and the vector order parameter to iden- 
tify the transitions. The average local-mode vector is 
proportional to the polarization. Note that, unless it is 
indicated, the MC simulations for the present effective 
Hamiltonian are performed at zero external pressure. 



IV. FINITE-TEMPERATURE RESULTS 



TABLE II: BaTiOa transition temperatures, in K, between 
cubic (C), tetragonal (T), orthorhombic (O), and rhombohe- 
dral (R) phases, as obtained from the shell model and from 
the effective Hamiltonian. The last three rows correspond to 
effective Hamiltonians modified as indicated in the text. In 
the first column, percentage error relative to the shell model 
is given in parentheses. 





C-T 


T-O 


O-R 


Shell model 


210 


135 


100 


Effective Hamiltonian 


150 (-28%) 


110 


85 


i/cff different a and 7 


150 (-28%) 


120 


100 


HcB + PcB 'by hand' 


185 (-12%) 


125 


95 


HcB + computed Peff 


165 (-21%) 


120 


90 



qucnce of transitions (cubic to tetragonal to orthorhom- 
bic to rhombohedral with decreasing T). However, the 
agreement is far from perfect for the Tc values, listed 
in Table |l[ Clearly the effective Hamiltonian underesti- 
mates the Tc's, especially for the C-T transition where 
the transition temperature is too low by ~30%. 

This shows that for an effective Hamiltonian of this 
form. Type II errors are quite significant, and in fact are 
comparable to the discrepancies found when comparing 
BaTiOa ab-initio effective-Hamiltonian results against 
experimental measurements on real BaTiOs samples. 
This strongly suggests that errors in first-principles cal- 
culations account for at most only part of the latter dis- 
crepancy. In the following, we will investigate the Type II 
errors in more detail and identify approaches to system- 
atic reduction of these errors, returning to the discussion 
of first-principles effective Hamiltonians in Sec. VII. 

V. ANALYSIS OF SOURCES OF THE 
DISCREPANCIES 

In this Section, we analyze three sources of error in 
the construction of the effective Hamiltonian that could 
possibly lead to the calculated underestimates of the Tc 
values. First, we focus on the specification of the ferro- 
electric mode unit vector, which determines the precise 
set of degrees of freedom described by the effective Hamil- 
tonian. Second, we consider the effect of the truncation 
of the Taylor expansion in the specified degrees of free- 
dom, with particular attention to the neglect of certain 
higher-order couplings within the effective Hamiltonian 
subspace. Third, we consider the consequences of omit- 
ting the higher-frequency phonon branches from the ef- 
fective Hamiltonian. 



Figure |l|a shows the three components of the mean po- 
larization as a function of temperature as obtained from 
the MC simulations of the effective Hamiltonian, while 
Fig. |l]b shows the corresponding results obtained directly 
from the shell-model MD simulations. It is apparent that 
the effective Hamiltonian correctly reproduces the se- 



A. Specification of ferroelectric local mode vector 

One of the first choices that was made in the construc- 
tion of the effective Hamiltonian was the detailed spec- 
ification of the ferroelectric local mode vector. As ex- 
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plained in Sec. Ill, we chose a Ti-centered displacement 
pattern selected in such a way that a uniform superpo- 
sition of these local mode displacements gives a periodic 
displacement pattern corresponding to the unstable (fer- 
roelectric) mode eigenvector of the force-constant matrix 
in the cubic structure. This is only one of many possible 
approaches, and questions may arise as to whether this 
is the best choice and how much difference it would make 
if we had made a different choice. 

One way of addressing these questions is to test how 
completely the chosen local mode vectors span the space 
of distortions that are actually encountered in the full 
atomistic finite-temperature shell-model simulation. We 
projected the shell-model MD trajectories at a given tem- 
perature onto the 12 optical zone-center normal modes 
(i.e., force-constant eigenvectors) of the cubic phase (four 
sets of three- fold degenerate modes). As expected, we 
found that the mode branches included in our effective 
Hamiltonian subspace account for almost all ('^90%) of 
the observed atomic displacements. This suggests that 
the approximation of keeping only these modes in the 
effective Hamiltonian is a good one. 

A second approach is to try a different procedure for 
defining the identity of the local mode vector. In par- 
ticular, one could think of a construction designed to 
optimize the description in the neighborhood of the fer- 
roelectric ground state. For instance, the local mode vec- 
tor could be fitted to the ground state of the system that 
is obtained when the atomic positions are fully relaxed. 
Such a procedure would effectively incorporate the ef- 
fect of the anharmonic couplings between included and 
excluded modes while not increasing the number of vari- 
ables considered in the model. In order to quantify the 
effect of this change, we assume that this alternative local 
mode definition mainly affects the anharmonic parame- 
ters in 

S^'^'f = K2U^ + au\ -f -i(u\^u\ -I- u\u\^ + u\^u\,). (2) 

We thus recalculated a and 7 exactly to reproduce the 
energy of the fully relaxed tetragonal energy minimum, 
and to get the best compromise for the energies of the 
orthorhombic and rhombohedral minima. By "fully re- 
laxed" we mean that the atomic positions were allowed to 
relax with the cell constrained to be the equilibrium cubic 
cell; this is consistent with the fact that we did not re- 
calculate any mode-strain coupling parameter. The new 
a and 7 are 0.811 and — 0.916 a. u. respectively, which 
are very similar to the values in Table |. The smallness 
of the correction reflects the fact that the fully relaxed 
energy minima are very close to those described by the 
original effective Hamiltonian, the differences being of 
the order of 0.01 mHa. Keeping all other parameters un- 
changed, we repeated the MC simulations at finite tem- 
perature and found that the transition teinperatures of 
the T-0 and 0-R transitions (see Table ||) are sensi- 
tive to these small changes in parameter values, giving 
a 10% improvement compared with our original effective 
Hamiltonian. However, the large discrepancy in the C-T 



transition temperature is unchanged. 

We therefore conclude that a change in the defini- 
tion of the local-mode displacement pattern is unlikely 
to be sufficient to eliminate the discrepancy between the 
effective-Hamiltonian and shell-model results. It is nec- 
essary, therefore, to look elsewhere. Nevertheless, the re- 
sults do show that the details of the fitting procedure can 
have a significant effect on the transition temperatures. 



B. Neglect of higher-order terms in the Taylor 
expansion 

We now return to our initial choice of relevant degrees 
of freedom, and ask whether the corresponding energy 
landscape is sufficiently well described by the truncated 
Taylor expansion that defines the effective Hamiltonian. 
The quadratic elastic energy E^^'^ is easily checked to be 
accurate. The dipole-dipole interactions in E'^^^ will be 
harmonic as long as the local polarization is linear in the 
atomic displacements, and this approximation is valid for 
BaTiOa. The terms that require further consideration 
are £;'^hort^ ^.int 

Higher-order terms in E^'^^'^ should aim at a better de- 
scription of the double- well potentials associated with the 
ferroelectric instabilities. We checked, however, that in- 
cluding sixth- and eighth-order terms does not improve 
the fit significantly. In particular, the well depths, which 
are the effective-Hamiltonian feature most directly re- 
lated to the value of the transition temperatures, are very 
well reproduced by the quartic A more accurate de- 

scription would yield energy wells around 1% shallower, 
which would probably lead to a very tiny decrease in the 
C-T transition temperature. 

Higher-order terms in iJ^^hort j-gpj-esent anharmonic 
couplings between neighboring local modes and would 
provide a correction to the local-anharmonicity approxi- 
mation. One can fit such terms by looking at the double- 
well potentials associated with the antiferroelectric insta- 
bilities of shell-model BaTiOs at zone-boundary points X 
and M. The fourth-order terms associated with such wells 
will be a combination of the parameters a and 7 in Eq. ^ 
and the new quartic parameters in E^"^™"^ . However, we 
find that these new quartic parameters are very small 
and can be safely neglected. More precisely, they consti- 
tute 0.05% and 5% of the total fourth-order term for the 
X and M instabilities, respectively, and result in slightly- 
deeper zone-boundary energy wells. Their probable ef- 
fect is a minor decrease in the transition temperatures, 
because of an enhanced competition between zone-center 
and zone-boundary instabilities. 

One could think of improving on E^"^™^ by including 
couplings between further neighbors (following Ref. |^, we 
included couplings up to third neighbors in our Hamilto- 
nian). This would allow us to improve the description of 
the dispersion branches of the relevant modes through- 
out the Brillouin zone. However, we checked that if our 
Hamiltonian is fitted including couplings only up to sec- 
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ond neighbors the transition temperatures change by less 
than 10 K. Hence, we can assume our Hamiltonian is weU 
converged in this respect. 

Finally, the description of the interaction between 
strain and local modes can be improved by including 
more terms in i?'"'. In particular, we have found that 
the r]iuf^ term is not negligible and would modify the ef- 
fective Hamiltonian so as to yield higher transition tem- 
peratures. Specifically, we find that the coefficient of the 
Vi'^ix term is negative and would thus favor a state in 
which the system polarizes and expands along a Carte- 
sian direction. However, in the next section we will see 
that the main Type H error has a different origin, and 
so we leave explicit consideration of this correction for 
future work. 



C. Effect of excluded modes: thermal expansion 

Finally, we consider the original decision to reduce the 
number of degrees of freedom in the effective Hamilto- 
nian to one vector degree of freedom per cell to repre- 
sent ferroelectric distortions. Even if an optimal set of 
local-mode variables is chosen (Sec. VA) and all neces- 
sary terms in the Taylor expansion are kept (Sec. VB), 
there still may be errors associated with this fundamen- 
tal approximation. For example, anharmonic couplings 
between included and excluded modes are neglected, as 
are anharmonic couplings between excluded modes and 
of excluded modes to strain. 

The effects of neglecting these anharmonic couplings 
are clearly seen in the thermal expansion. In fact, in 
raising the temperature from to 300 K in our simula- 
tions, we find that the volume of the shell- model system 
increaseshy 0.4%, while that of the cffcctive-Hamiltonian 
system decreases by 0.6%. This indicates that the effec- 
tive Hamiltonian treatment of the thermal expansion is 
qualitatively incorrect. Moreover, given the well-known 
sensitivities of the transition temperatures to volume, 
this effect could be substantial. Moreover, it correctly 
predicts that we would underestimate transition tem- 
peratures, since they are reduced at smaller lattice con- 
stants. 

To check whether the thermal expansion effect is re- 
sponsible for the dominant errors in Tc, we made the 
following test. We completely eliminated the volume ef- 
fect by carrying out both simulations at a fixed volume 
of (4.012 A)"^. Using the shell model, we have found Tc 
values of 190, 130, and 95 K, while the corresponding 
effective-Hamiltonian values are 180, 125, and 100 K, re- 
spectively. The error in the C-T transition temperature, 
which was around 60 K in the zero-pressure simulations, 
is reduced to ^10 K. 



D. Summary 

We thus arrive at the important conclusion that the 
poor description of thermal expansion effects is the dom- 
inant source of error in the effective-Hamiltonian descrip- 
tion. These shortcomings in the description of thermal 
expansion, and some preliminary attempts to correct for 
them, will be described in the following Section. Smaller 
errors (probably amounting to no more than 5-10% of the 
Tc values) are associated with the other approximations 
discussed in Subsecs. VA-B. 



VI. IMPROVED TREATMENT OF THERMAL 
EXPANSION 

Given the conclusion of Sec. V.D, we are strongly moti- 
vated to improve the effective-Hamiltonian treatment of 
thermal expansion. First, we investigate the thermal ex- 
pansion in more detail. Figure ^ shows the pseudocubic 
lattice parameter a = V^^^ as a function of tempera- 
ture as predicted by the shell model (full circles) and by 
the effective Hamiltonian (full triangles). Both models 
exhibit volume anomalies ("kinks") at the ferroelectric 
phase transition temperatures, indicative of their first- 
order character. However, the overall trends in volume 
vs. temperature are quite different. 

The thermal expansion displayed by the shell model, 
after a proper rescaling of temperature and_cell parame- 
ter, closely resembles that of real BaTiOa.EEl This virtue 
of the model is related to the fact that it includes all 
the degrees of freedom of the system and a sufficiently 
accurate description of their relevant anharmonicities. 
The effective Hamiltonian, on the other hand, does not 
properly account for the thermal expansion of the sys- 
tem, and actually leads to a contraction with increasing 
temperature in the range of the polar phases. The rea- 
son for such a contraction is that the volume is strongly 
coupled to the magnitude of the local dipoles and, as 
these tend to decrease with increasing temperature as 
the paraelectric phase is approached, the volume tends 
to decrease as well. Equivalently, the thermal contrac- 
tion can be attributed to negative Griineisen parameters 
associated with portions of the relevant branches; these 
are overwhelmed by positive contributions from higher 
modes in the shell-model system, but not in the effective- 
Hamiltonian system where the higher modes are absent. 

We next ask what happens if the effective-Hamiltonian 
simulations are carried out with a cell volume that is 
constrained 'by hand' to have the correct dependence 
on temperature as given by the shell-model system. A 
simple way of doing this in practice is to apply a (neg- 
ative) external pressure pcff to the effective-Hamiltonian 
system; this fictitious pressure can be thought of as aris- 
ing from the thermal expansion effects of the excluded 
modes. We implement this approximately by taking 
PcB to be linear in temperature in such a way that the 
effective-Hamiltonian equilibrium volume coincides with 



7 



4.015 



g 4.01 
"o 

> 

4.005 



Shell model 



2v\. A A A- 



p 'by hand' 



computed 



Effective Hamiltonian 



the classical shell model, PcS takes the form 



50 100 150 200 250 300 

Temperature (K) 



FIG. 2; Pseudocubic lattice parameter a = V^^^ {V = cell 
volume) for BaTiOs as a function of the temperature as pre- 
dicted by the effective Hamiltonian (full triangles) and by the 
shell model (full circles). Open triangles and squares corre- 
spond to effective Hamiltonian results under an external pres- 
sure adjusted 'by hand' and ab-imtio respectively (see text for 
details) . 



the shell-model one at two temperatures, taken to be 
10 and 300 K, bracketing the relevant range. We then 
find that pcff is required to be — l.SGPa at 300 K while 
almost vanishing at 10 K. The results of MC simulations 
under this external pressure are presented in Fig. ^ (open 
squares); the values of the transition temperatures are 
listed in Table ||. The improvement in the agreement 
with the shell-model calculations, in particular in the case 
of the C-T transition temperature, is remarkable. 

However, such an ad-hoc approach is not consistent 
with the spirit of first-principles based approaches; one 
would prefer a way to calculate the effective pressure PcS 
ab-initio. We have attempted to do so by employing the 
so-called quasiharmonic approximation (see Chapter 25 
of Ref. ^l|) . Within this approximation, the pressure that 
develops in a harmonic crystal with volume-dependent 
phonon frequencies is (see Eq. (25.5) of Ref. El]) 




(3) 



where is the equilibrium energy of the system, lOs (k) 
is the phonon frequency of branch s at point k of the 
Brillouin zone (BZ), and the summations run over all 
branches and k points. Now the pressure exerted by the 
excluded modes is obtained from Eq. (^ by removing 
U'^'^ and restricting the sums to the excluded modes s'. 
Taking the classical limit ?i ^ in order to compare with 



PeS 



(4) 



which is linear in temperature and proportional to the 
volume derivative of the phonon frequencies. 

We must be cautious about the approximations in- 
volved in the use of Eq. (^ or its quantum-mechanical 
version. The quasiharmonic approach is not really well 
suited to deal with phase transitions, which are strongly 
anharmonic phenomena. Using it in the present con- 
text relies on the assumption that the excluded modes 
(more precisely, the volume derivatives of their frequen- 
cies) are not significantly affected by the strong fluctua- 
tions and phase transitions associated with the relevant 
local modes. 

In order to assess the utility of the quasiharmonic ap- 
proach here, we have focused on the cubic-to-tetragonal 
transition and calculated Poff using the volume depen- 
dences of the exckijded-mode frequencies in the cubic 
paraelectric phase.EHl We find that the BZ sum in Eq. (Q) 
can be evaluated with good accuracy using information 
from the high-symmetry k-points only, and that the sum 
of logarithms of the hard-mode frequencies depends lin- 
early with volume in the relevant volume range, thus al- 
lowing us to take pcs as independent of volume. The 
Pcff calculated in this way shows improved agreement 
with the exact shell-model results as regards both the 
transition temperatures (denoted by "computed Pcff" in 
Table |l|) and the thermal expansion (open triangles in 
Fig. ^). However, the results are still far from satisfac- 
tory, with a substantial error remaining for the C-T tran- 
sition temperature. These discrepancies are probably 
connected with the shortcomings of the quasiharmonic 
approximation discussed above. Further investigations 
along these lines are clearly needed, but fall beyond the 
scope of the present paper. 

In summary, the proposed correction based on the 
quasiharmonic approximation of Eq. (^) accounts cor- 
rectly for only a fraction (perhaps a third) of the thermal- 
expansion error. Unfortunately, then, we are not yet in 
a position to propose a fully ab initio approach to the 
thermal expansion problem in the context of effective- 
Hamiltonian methods. 



VII. CONCLUSIONS 

The main weakness of our effective-Hamiltonian de- 
scription of shell- model BaTiOa is the poor description of 
the thermal expansion. Focusing on the cubic to tetrago- 
nal transition temperature (Tct), we have found that the 
effective Hamiltonian produces a 28% error, while we can 
estimate from our constant- volume calculations that this 
error should be reduced to 5% if the thermal expansion 
were properly modeled. We have also seen that including 
the thermal expansion 'by hand' allows us to reduce the 
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error to about 12%; in other words, this correction seems 
to account for 70% of the total error associated with the 
thermal expansion. 

It is tempting to apply these same percentages in order 
to estimate the sources of error arising in the compari- 
son of the first-principles effective Hamiltonian transi- 
tion temperatures with real experiment. However, this 
should be done cautiously. For example, anharmonicities 
or thermal-expansion effects might either be exaggerated 
or underestimated by the shell model. With this in mind, 
we consider the effective-Hamiltonian study of BaTiOs 
by Zhong et al. which led to Tct = 300 K, 25% be- 
low the experimental value of 400 K. This was a classical 
calculation; should quantum effects be considered, the 
calculated Tct would be smaller by about 30Kjl3 and 
thus the error would go up to ^^30%. 

We performed classical MC simulations with the effec- 
tive Hamiltonian of Zhong et al. including the thermal 
expansion of the system 'by hand' under the condition 
that the computed volume should coincide with the ex- 
perimental one at T = 473 K. This resulted in an error 
of 18% in Tct, which would become ~25% if we include 
the estimated quantum effects. Hence, it seems that the 
improvement in this case is not as large as it was for the 
effective Hamiltonian fitted to shell- model BaTiOs. If 
we follow what we have learned from the case of shell- 
model BaTiOa and assume that including the thermal 
expansion 'by hand' corrects 70% of the total thermal- 
expansion error, we can estimate that a quantum first- 
principles effective-Hamiltonian calculation with perfect 
thermal expansion would still result in a 20% underesti- 
mate of Tct- It seems reasonable to assume that Type II 
errors other than the thermal expansion, as well as de- 
tails of the fitting procedure, are responsible for a further 
5% error in Tct- This suggests that a calculation free of 
Type II and Type IV errors would yield a Tct that would 
still be about 15% below experiment. While this line of 



reasoning is tenuous, we nevertheless believe it gives the 
best current estimate for the magnitude of the error that 
should be attributed to the first-principles methods used 
to construct the effective Hamiltonian (in particular, the 
LDA). 

In summary, in this work we have analyzed the errors 
associated with the first-principles effective Hamiltonian 
method that has been developed for the treatment of 
the thermodynamics of perovskite ferroelectrics. More 
specifically, by considering the effective-Hamiltonian de- 
scription of the shell-model for BaTiOa of Tinte et ai, 
we have been able to isolate and study in detail the er- 
rors intrinsic to the effective-Hamiltonian approximation 
(Type II errors). We have found that the main Type II 
error is associated with a poor description of the thermal 
expansion of the system. We have discussed an easily- 
implemented first-principles correction that takes into 
account some contributions of excluded modes. Unfor- 
tunately, this scheme seems to account for only about a 
third of the total thermal-expansion error. More elab- 
orate schemes (involving a more thorough treatment of 
the couplings of the phonon modes to each other and to 
strain and polarization) might substantially reduce the 
error, but it remains for the future to explore and imple- 
ment such schemes. Finally, we have argued that in the 
case of the comparison of the first-principles effective- 
Hamiltonian calculations on BaTiOs with real experi- 
ment. Type II errors do not seem to be responsible for the 
entire discrepancy. Our results suggest that the Type I 
errors associated with the use of the LDA and other first- 
principles technicalities may be of the same magnitude as 
the thermal-expansion error. 
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